Proyecto 1 computer vision

Miguel López

ID 1001014378

Image description

Convocatoria 1 - Proyecto 1¶

0) Cargar una de las imágenes histológicas¶

In [33]:
# Utilizar la librería skimage.io para leer la imagen 'histo_x.jpg' en formato RGB.
# Normalizar la imagen para que los píxeles se encuentren en el rango [0, 1]
# Visualizar la imagen
import skimage.io as io
import matplotlib.pyplot as plt

image = io.imread('imgs/histo_1.jpg')

print("Valor máximo antes de normalizar",image.max())

image = image / 255

print("Valor máximo después de normalizar",image.max())

print(image.shape)

plt.imshow(image)
plt.show()
Valor máximo antes de normalizar 255
Valor máximo después de normalizar 1.0
(1024, 1024, 3)
No description has been provided for this image

1) Realizar una transformación de color para convertir la imagen al espacio de color CMYK¶

In [34]:
# Extraer la componente magenta de la imagen (que corresponde a la región tisular)
# Visualizar la imagen del canal magenta

#Al ser RGB, el magenta es la combinación de rojo, azul y un poco de verde

rojo = image[:,:,0]
verde = image[:,:,1]
azul = image[:,:,2]

mascara = (rojo > 0.5) & (verde < 0.5) & (azul > 0.5)

magenta = image.copy()

magenta[~mascara] = 0

plt.title("Magenta")
plt.imshow(magenta, cmap='gray')
plt.show()
No description has been provided for this image
In [35]:
def pink_extraction(image):

    rojo = image[:,:,0]
    verde = image[:,:,1]
    azul = image[:,:,2]

    mascara = (rojo > 0.5) & (verde < 0.5) & (azul > 0.5)

    magenta = image.copy()

    magenta[~mascara] = 0

    plt.title("Magenta")
    plt.imshow(magenta, cmap='gray')
    plt.show()
    return magenta

magenta = pink_extraction(image)
No description has been provided for this image
In [36]:
#Extracción directa de treshold
import cv2

image_uint8 = (image * 255).astype('uint8')

img = cv2.cvtColor(image_uint8, cv2.COLOR_RGB2GRAY)

plt.title("Histograma")
plt.hist(img.flatten(), bins=256, range=[0,256])
plt.show()

plt.title("Imgagen en escala de grises")
plt.imshow(img, cmap='gray')
plt.show()

mascara_magenta = (img > 200).astype('uint8')
mascara_magenta = ~mascara_magenta

plt.title("Máscara")
plt.imshow(mascara_magenta, cmap='gray')
plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

2) Umbralizar la imagen para separar los píxeles del fondo de la región tisular¶

In [37]:
# Aplicar un filtro gaussiano de tamaño 5x5 y después utilizar el método de Otsu de manera que
# los píxeles correspondientes al lumen y al background de la imagen sean 1s y el resto de los píxeles tengan un valor de 0.
# Nota: Recordar que el método de Otsu requiere como input una imagen en el rango [0-255] en formato "uint8".
# Visualizar la máscara resultante

from cv2 import GaussianBlur

def zone_extraction(img_path):
    image = io.imread(img_path)

    plt.title("Imagen original")
    plt.imshow(image)
    plt.show()

    image_uint8 = (image).astype('uint8')

    img = cv2.cvtColor(image_uint8, cv2.COLOR_RGB2GRAY)

    gauss = GaussianBlur(img, (5, 5), 0)

    plt.title('Filtro Gaussiano')
    plt.imshow(gauss, cmap='gray')
    plt.show()

    t, num_mask = cv2.threshold(gauss, 0, 1,cv2.THRESH_OTSU)
    #Se convierte a boleano
    mask = num_mask == 1

    inv_mask = (~mask * 255).astype("uint8")


    plt.title('Máscara Otsu t=' + str(t))
    plt.imshow(mask, cmap='gray')
    plt.show()

    filter_img = cv2.bitwise_and(image, image, mask=inv_mask)
    plt.title("Imagen original con filtro inverso")
    plt.imshow(filter_img)
    plt.show()

    return mask, img, image

mask, img, original_image = zone_extraction("imgs/histo_1.jpg")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [ ]:
 

3) Limpiar la imagen eliminando los artefactos de lumen (objetos blancos pequeños que no son lúmenes)¶

In [38]:
# Utilizar la librería skimage.morphology.remove_small_objects para eliminar aquellos objetos cuya área sea menor a 300 píxeles
# Más información en https://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.remove_small_objects
# Visualizaer la máscara resultante

from skimage import morphology


mask_small = morphology.remove_small_objects(mask, 300)

#Aplicamos la máscara a la imagen
plt.title("Eliminación de objetos pequeños")
plt.imshow(mask_small, cmap="gray")
plt.show()
No description has been provided for this image

4) Rellenar con 0s el fondo de la imagen para quedarnos únicamente con los lúmenes¶

In [39]:
# Aplicar el algoritmo de expansión a partir de semillas (region growing) de manera que únicamente los lúmenes sean blancos
# y el resto de la imagen negra. Pista: utilizar dos semillas. Nota: Se pueden fijar las semillas de manera manual, pero
# se valorará positivamente a aquell@s que desarrollen una función para encontrarlas automáticamente.
# Visualizar la máscara resultante.
from skimage import measure
import numpy as np

def find_seeds(mask, min_area):

    labeled = measure.label(mask * 1)
    labeled = cv2.convertScaleAbs(labeled)

    seeds = []
    contours, _ = cv2.findContours(labeled, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

    #De esta forma obtenemos los centroides de las regiones identificadas.
    areas = []
    for contour in contours:
        area = cv2.contourArea(contour)
        areas.append(area)
        #Luego de examinar las areas se agrega un filtro para solo colocar
        #semillas en los lugares deseados
        if area >= min_area:
            border_pixels = contour.squeeze()
            if len(border_pixels) > 0:
                seed = tuple(border_pixels[0])
                seeds.append(seed)
                continue


    print(sorted(areas, reverse=True))
    return seeds, labeled

def apply_flood_fill(img, seeds):

    h, w = img.shape
    ref = np.zeros((h+2, w+2), np.uint8)

    img_copy = img.copy()

    img_copy = img_copy.astype("uint8")

    # seeds = [(0, 0),(1000, 1000)]

    for seed in seeds:
        cv2.floodFill(img_copy, None, seed, 0)

    # cv2.floodFill(img_copy, None, (0, 0), 255)

    img_copy = cv2.bitwise_not(img_copy)

    return img_copy

seeds, labeled = find_seeds(mask_small, 11e4)

new_mask = apply_flood_fill(mask_small * 255, seeds)
new_mask = cv2.bitwise_not(new_mask)

img_copy = img.copy()

for seed in seeds:
    cv2.circle(img_copy, tuple(seed) , 5, 255, thickness=20)

print("Semillas encontradas en", seeds)
plt.title("Semillas encontradas")
plt.imshow(img_copy, cmap="gray")
plt.show()

plt.title("Obtención de solo los lumens")
plt.imshow(new_mask, cmap="gray")
plt.show()

plt.title("Subsecciones de la mascara")
plt.imshow(labeled, cmap="gray")
plt.show()
[316521.5, 115006.5, 12884.0, 12008.0, 6076.5, 3305.0, 2161.5, 1710.0, 1263.0, 1152.0, 1000.0, 767.5, 744.0, 590.0, 383.0, 325.5, 318.0]
Semillas encontradas en [(np.int32(975), np.int32(615)), (np.int32(0), np.int32(0))]
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

5) Rellenar los objetos de los lúmenes¶

In [40]:
# Rellenar los lúmenes con la función binary_fill_holes de la librería scipy.ndimage.morphology
# Visualizar la máscara resultante
import scipy as sp

mask_no_holes = sp.ndimage.binary_fill_holes(new_mask)

plt.title("Subsecciones de la mascara")
plt.imshow(mask_no_holes, cmap="gray")
plt.show()
No description has been provided for this image
In [41]:
def lumen_separation(mask, min_area, img):

    mask_small = morphology.remove_small_objects(mask, 300)

    seeds, labeled = find_seeds(mask_small, min_area)

    new_mask = apply_flood_fill(mask_small *255, seeds)
    new_mask = cv2.bitwise_not(new_mask)

    img_copy = img.copy()

    for seed in seeds:
        cv2.circle(img_copy, tuple(seed) , 5, 125, thickness=20)

    print("Semillas encontradas en", seeds)
    plt.title("Semillas encontradas")
    plt.imshow(img_copy, cmap="gray")
    plt.show()

    mask_no_holes = sp.ndimage.binary_fill_holes(new_mask)


    plt.title("Eliminación de objetos pequeños")
    plt.imshow(mask_small, cmap="gray")
    plt.show()

    plt.title("Obtención de solo los lumens")
    plt.imshow(new_mask, cmap="gray")
    plt.show()

    plt.title("Subsecciones de la mascara")
    plt.imshow(mask_no_holes, cmap="gray")
    plt.show()

    return  mask_no_holes

mask_no_holes = lumen_separation(mask, 11e4, img)
[316521.5, 115006.5, 12884.0, 12008.0, 6076.5, 3305.0, 2161.5, 1710.0, 1263.0, 1152.0, 1000.0, 767.5, 744.0, 590.0, 383.0, 325.5, 318.0]
Semillas encontradas en [(np.int32(975), np.int32(615)), (np.int32(0), np.int32(0))]
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

6) Detectar y dibujar los contornos de los lúmenes sobre la imagen original¶

In [42]:
# Dibujar los contornos de los lúmenes en color verde sobre la imagen original RGB. Nota: Utilizar los flags necesarios
# para que los contornos en verde sean perfectamente visibles. 
# Visualizar la imagen superpuesta

mask_no_holes_u = (mask_no_holes * 255).astype("uint8")


canny = cv2.Canny(mask_no_holes_u, 255/3, 255)
kernel = np.ones((2,2), 'uint8')
canny = cv2.dilate(canny, kernel, iterations=3)

#Mejoramos la calidad de los bordes


contours, _ = cv2.findContours(canny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)

image_with_contours = cv2.drawContours(image_uint8.copy(), contours, -1, (0, 255, 0), 2)

plt.imshow(image_with_contours)
plt.title('Image with Contours')
plt.show()
No description has been provided for this image

7) Identificar y cropear el lumen más grande¶

In [43]:
# Determinar cuál es el lumen de mayor área y hacer un crop del mismo sobre la imagen original RGB.
# Visualizar el lumen cropeado.

#Se identifica el lumen más grande
largest_contour = max(contours, key=cv2.contourArea)

#Se realiza un crop rectangular de este para tener un contexto de la imagen

x, y, w, h = cv2.boundingRect(largest_contour)

cropped_lumen = image[y:y+h, x:x+w]
cropped_mask = mask_no_holes_u[y:y+h, x:x+w]

plt.imshow(cropped_lumen)
plt.title("Mayor lumen")
plt.show()
No description has been provided for this image
In [44]:
def big_lumen(mask_no_holes, img):
    mask_no_holes_u = (mask_no_holes * 255).astype("uint8")


    canny = cv2.Canny(mask_no_holes_u, 255/3, 255)
    kernel = np.ones((2,2), 'uint8')
    canny = cv2.dilate(canny, kernel, iterations=3)

    #Mejoramos la calidad de los bordes

    contours, _ = cv2.findContours(canny, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)


    image_with_contours = cv2.drawContours(img, contours, -1, (0, 255, 0), 2)

    plt.imshow(image_with_contours)
    plt.title('Contornos en la imagen')
    plt.show()

    #Se identifica el lumen más grande
    largest_contour = max(contours, key=cv2.contourArea)

    #Se realiza un crop rectangular de este para tener un contexto de la imagen

    x, y, w, h = cv2.boundingRect(largest_contour)

    cropped_lumen = img[y:y+h, x:x+w]
    cropped_mask = mask_no_holes_u[y:y+h, x:x+w]

    plt.imshow(cropped_lumen)
    plt.title("Mayor lumen")
    plt.show()

    return cropped_mask, cropped_lumen


cropped_mask, cropped_lumen = big_lumen(mask_no_holes, original_image)
No description has been provided for this image
No description has been provided for this image

8) Extraer 13 características geométricas que permitan caracterizar el lumen recortado¶

In [45]:
# Calcular las siguientes características del crop del lumen de mayor área, redondeando su valor hasta el cuarto decimal.
# 1) Área
# 2) Área de la bounding box
# 3) Área convexa
# 4) Exentricidad
# 5) Diámetro equivalente
# 6) Extensión
# 7) Diámetro Feret
# 8) Longitud del eje mayor
# 9) Longitud del eje menor
# 10) Orientación
# 11) Perímetro
# 12) Solidez
# 13) Compacidad


area = cv2.contourArea(largest_contour)

# Calculate bounding box area
bounding_box_area = w * h

# Calculate convex hull and convex area
hull = cv2.convexHull(largest_contour)
convex_area = cv2.contourArea(hull)

# Calculate eccentricity, major axis length, minor axis length, orientation, and solidity using regionprops
labeled_mask = measure.label(cropped_mask)
props = measure.regionprops(labeled_mask)

# Find the properties of the largest region
largest_region = max(props, key=lambda r: r.area)

eccentricity = largest_region.eccentricity
major_axis_length = largest_region.major_axis_length
minor_axis_length = largest_region.minor_axis_length
orientation = largest_region.orientation
solidity = largest_region.solidity

# Calculate equivalent diameter
equivalent_diameter = np.sqrt(4 * area / np.pi)

# Calculate extent
extent = area / bounding_box_area

# Calculate Feret diameter (maximum distance between any two points on the contour)
feret_diameter = max([np.linalg.norm(p1 - p2) for p1 in largest_contour for p2 in largest_contour])

# Calculate perimeter
perimeter = cv2.arcLength(largest_contour, True)

# Calculate compactness
compactness = (perimeter ** 2) / (4 * np.pi * area)

# Print the results rounded to the fourth decimal place
print(f"Area: {area:.4f}")
print(f"Bounding Box Area: {bounding_box_area:.4f}")
print(f"Convex Area: {convex_area:.4f}")
print(f"Eccentricity: {eccentricity:.4f}")
print(f"Equivalent Diameter: {equivalent_diameter:.4f}")
print(f"Extent: {extent:.4f}")
print(f"Feret Diameter: {feret_diameter:.4f}")
print(f"Major Axis Length: {major_axis_length:.4f}")
print(f"Minor Axis Length: {minor_axis_length:.4f}")
print(f"Orientation: {orientation:.4f}")
print(f"Perimeter: {perimeter:.4f}")
print(f"Solidity: {solidity:.4f}")
print(f"Compactness: {compactness:.4f}")
Area: 15535.0000
Bounding Box Area: 39105.0000
Convex Area: 28602.0000
Eccentricity: 0.8445
Equivalent Diameter: 140.6406
Extent: 0.3973
Feret Diameter: 246.4650
Major Axis Length: 231.6881
Minor Axis Length: 124.0865
Orientation: 0.9970
Perimeter: 1117.6022
Solidity: 0.4872
Compactness: 6.3981
In [46]:
import math
import pandas as pd


def get_metric_data(cropped_mask):

    new_lab, new_num =measure.label(cropped_mask, return_num=True)


    # Crear un DataFrame para almacenar las características geométricas
    columns = ['ID', 'Area', 'BBox Area', 'Convex Area', 'Eccentricity', 'Equiv Diameter', 'Extent', 'Major Axis', 'Minor Axis', 'Orientation', 'Perimeter', 'Solidity', 'Compactness', 'Rectangularity']
    data = []

    for i in range(1, new_num+1):
        objeto = new_lab == i
        prop = measure.regionprops(objeto.astype(np.uint8))[0]

        features = {
            'ID': i,
            'Area': np.round(prop.area, 4),
            'BBox Area': np.round(prop.bbox_area, 4),
            'Convex Area': np.round(prop.convex_area, 4),
            'Eccentricity': np.round(prop.eccentricity, 4),
            'Equiv Diameter': np.round(prop.equivalent_diameter, 4),
            'Extent': np.round(prop.extent, 4),
            'Major Axis': np.round(prop.major_axis_length, 4),
            'Minor Axis': np.round(prop.minor_axis_length, 4),
            'Orientation': np.round(prop.orientation, 4),
            'Perimeter': np.round(prop.perimeter, 4),
            'Solidity': np.round(prop.solidity, 4),
            'Compactness': np.round(4 * math.pi * prop.area / prop.perimeter**2, 4),
            'Rectangularity': np.round(prop.area / prop.bbox_area, 4)
        }
        data.append(features)

    df = pd.DataFrame(data, columns=columns)
    return df

get_metric_data(cropped_mask)
Out[46]:
ID Area BBox Area Convex Area Eccentricity Equiv Diameter Extent Major Axis Minor Axis Orientation Perimeter Solidity Compactness Rectangularity
0 1 13379.0 37674.0 27463.0 0.8445 130.5169 0.3551 231.6881 124.0865 0.997 1161.9697 0.4872 0.1245 0.3551

2. Pipeline para hist.2¶

In [47]:
def pipeline(img_path, min_area):
    mask, img, original_image = zone_extraction(img_path)
    magenta = pink_extraction(original_image)
    mask_no_holes = lumen_separation(mask, min_area, img)
    cropped_mask, cropped_lumen = big_lumen(mask_no_holes, original_image)
    df = get_metric_data(cropped_mask)

    return df, cropped_mask, cropped_lumen, img


df2, cropped_mask2, cropped_lumen2, img2 = pipeline("imgs/histo_2.jpg", 3e3)
df2
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
[843344.5, 2194.5, 1966.0, 974.0, 884.0, 874.5, 730.0, 684.0, 604.0, 587.5, 483.5, 395.0, 367.5, 342.5, 330.0]
Semillas encontradas en [(np.int32(573), np.int32(0))]
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
Out[47]:
ID Area BBox Area Convex Area Eccentricity Equiv Diameter Extent Major Axis Minor Axis Orientation Perimeter Solidity Compactness Rectangularity
0 1 4675.0 7392.0 4875.0 0.8392 77.1518 0.6324 104.9016 57.0514 -0.5622 287.5635 0.9590 0.7104 0.6324
1 2 385.0 578.0 399.0 0.8947 22.1404 0.6661 34.8248 15.5540 -0.2476 88.6274 0.9649 0.6159 0.6661